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We present a formalism and numerical results for the energy loss of a charged particle scattered 
c I ^ . at an arbitrary angle from epitaxially grown multilayer graphene (MLG). It is compared with that 

' of free-standing graphene layers. Specifically, we investigated the effect of the substrate induced 

^.^1 energy gap on one of the layers. The gap yields collective plasma oscillations whose characteristics 

, are qualitatively and quantitatively different from those produced by Dirac fermions in gapless 

graphene. The range of wave numbers for undamped self-sustaining plasmons is increased as the 
gap is increased, thereby increasing and red-shifting the MLG stopping power for some range of 
, charged particle velocity. We also applied our formalism to interpret several distinct features of 

experimentally obtained electron energy loss spectroscopy (EELS) data. 

CJ ■ I. INTRODUCTION 

Electron energy loss spectroscopy (EELS) has been considered for many systems dating back to the classic paper of 
y Ritchie ^ on surface plasmons for a slab of dielectric material in the local limit and subsequently generalized to a non- 
^ local dielectric function by Gumbs and Horing.— The non-locality was included in Ref.[3 through the screening of the 
^ , electron-electron interaction produced by single-particle excitations and plasmon modes. Fessatidis, et al.— recently 
^ ■ adopted the formalism of Horing, et al.— for a two-dimensional electron gas (2DEG) and Gumbs and Balassis— for a 
^ ' nanotube to graphene with the aid of the polarization function calculated by Wunsch, et al. ^ for conventional Dirac 
electrons. 

There have been several recent papers which considered the effects which a circularly polarized electromagnetic field 
Q (CPEF), spin-orbit interaction (SOI) ^ in suspended graphene or the sublattice symmetry breaking (SSB) by an 
underlying polar substrate—^— in epitaxial graphene may have on the energy band structure and plasma excitations— 
of a graphene sheet. Under these conditions, there is a gap between the valence and conduction bands as well as 
I ' between the intra-band and inter-band electron-hole continuum of the otherwise semimetallic Dirac system.— li^ii^ 
^ \ Additionally, the interplay between the single-particle excitations in the long wavelength limit results in dielectric 
. screening of the Coulomb interaction which produces an undamped plasmon mode that appears in the gap separating 
' the two types of electron-hole modes forming a continuum. By this we mean that a self-sustained collective plasma 
OO mode is not supported by exciting either the intra-band or inter-band single-particle modes only. Furthermore, since 
' the Dirac electrons now acquire a non-zero effective mass, one can produce a long wavelength plasmon mode as in the 
2DEG by intra-band excitations only. These properties of the plasma modes give rise to noticeable differences in the 
behavior of the stopping power of gapped graphene compared with conventional free-standing exfoliated multilayer 
\ graphene. 

" , [ • Multi-layer epitaxially grown graphene (MLG) may become a valuable and relatively cheap alternative to rather ex- 
J> I pensive exfoliated graphene. Recent angle-resolved photoemission spectroscopy (ARPES) experiments unambiguously 
demonstrated almost perfect Dirac cones on most of the layers. That is the electrons in the layers behave as if the 
layers are uncoupled in contrast to Bernal stacking in bi-layer graphene. Although some layers may establish bi-layer 
i structures the interaction between the first and the buffer layer (the one sitting on top of the SiC substrate) is always 
weak.— lii The spectral function of MLG on SiC substrate suggests the energy gap of few hundreds of meV. However 
the exact gap opening mechanism in epitaxial graphene is still under debate.— The complexity comes from the fact 
that the graphene sample sits on top of a buffer layer which provides an additional mid gap level, thus obscuring the 
exact energy dispersion curve and requires numerical ab initio calculations.— In addition, the DOS around Fermi 
energy failed to indicate the gap.— This ambiguity stimulated discussion regarding the symmetry breaking gap— 
versus the effects due to electron-electron interaction. ^1 This also indicates the importance of alternative techniques 
capable of identifying the induced gap by a buffer layer on a substrate. 

EELS may be employed to ascertain the plasmon frequencies in single and double layer graphene. The Raman 
shift of the scattered electrons provides both particle-hole and plasmon excitation frequencies, which are usually 
characterized by their spectral weight, a quantity that depends on the transferred energy huj and in- plane momentum 
hq. This allows mapping of the plasmon dispersions a;p(g).^^ Here we calculate EELS spectra in order to pinpoint 
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the effect of the substrate induced gap in MLG. 

The outhne of the remainder of this paper is as follows. In Sec.|lTl we derive the energy loss for charged particles 
incident on multiple electron layers at an arbitrary angle of incidence. In Sec. lIIIi this formalism is applied to MLG 
with a gapped buffer layer. Our numerical simulations are compared with experimental results of Lu and Loh^^ in 
Sec.|lVl A summary is presented in Sec. El 



II. ENERGY LOSS FORMALISM 



We shall consider an inhomogeneous medium (in the z direction) described by a nonlocal dielectric function e (1; 2). 
Here, we have introduced space-time points, such as i = {zi,ri,ti), where the z spacial coordinate has been separated 
out due to the inhomogeneity and making cylindrical coordinates suitable for our discussion. Our focus is the 
calculation of energy loss of a charged particle Ze moving along a prescribed trajectory zi + vti. Due to non- locality 
of the dielectric function, the charged particle experiences a frictional force given by the gradient of the effective 
potential as 

F(l) = -Ze[ViCi>tot(l)]i=,,+.t, . (1) 
The total potential can be determined from the following system of equations 



$tot(l) = Jd^2 e-i(l;2)$ext(2) , (2) 

V?$ext(l) = -— <^(1-21-Vti) . (3) 

Solving Eq. ([3]) by Fourier transformation and substituting the result into Eq. ^ and then into ([T}, we obtain the 
components of the frictional force as 



2 oo oo 

F||(l) / dz2 /^q||d'q|| / dq, e'^^^'-^'^^'^^'^W^e-^ {zi + v^h, Z2;qii,<l- v) , (4) 

(27r) €0 J J J 

— oo — oo 

2 oo oo 

^^^^^^"(^^ / I / ^g-"'"'^"'"'"*'^ '^'^ d (zi I v^h) + ^^^1^ q ■ v) . (5) 

— 00 —00 

We assume that the particle trajectory begins at i = — oo and ends at i = oo with its motion in vacuum. The net 
energy lost due to the frictional force between it and the plasma during this motion may be expressed as 



00 

W = J dtiv(ti)-F(l) . (6) 

— 00 

It is convenient to partition the loss a,s W — + W±. The energy loss due to parallel and perpendicular motions 
may be written as 



2 00 00 



M/i|(zi,v) = --^^^ / dh I dz2 I d\\\ (qil-vii) / dg. e*«^(^^--i-''^*i)g^2^-i(zi+Mi,Z2;g||,q-v) , (7) 



(27rreo 
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In Eq. (IS]), integration by parts was carried out, along with the condition that the charged particle starts and ends in 
vacuum: (—00, Z2; (Jn , q • v) — (00, Z2; , q • v) — 1. The above general expressions can be simplified further, 
for the practically important cases of parallel, perpendicular and traversing motion. 
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A. Charged Particle moving parallel to the surface (v^ = 0) 

First, we assume that the charged particle travels parallel to the x, y plane with velocity coraponent Vz — 0. In this 
case, to simplify Eq. ([7]), we can use the following identities: 

dq, 9-2e^9.(.2-.i) = 2n , (9) 

2911 

— C30 

OO 1 

, \ ,2 -W \ r, /" 7 2 /" COS0 dfcOS^) / /,\ / r^\ 

(q|| • V||j d q||e (zi, Z2; g|| , q|| ■ V|| j = -2z / / e [zi,Z2;q\\,q\\v\\cosO) (10) 



-1 



Vl — COS^ y 

^ - J (fqH (q|| • V||) 3me"^ (zi,Z2;g||,q|| • V||) . 
In the above Eq. pUj) . to integrate over the angle we used the parity of the inverse dielectric function, i.e., 

e"^ (zi,Z2;-q||,-a;) = [e^i (zi,Z2;q||,w)]* . (11) 
Consequently, we obtain the energy loss in the parallel case as 



„ 00 00 

W|| (zi,W||) = j y dti y dz2 y d q|| (qr^n) 2eog|| — ^i^^" (^i'^2;9||,q|| -vn) . (12) 
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Since the integrand in the above expression does not depend on time, it is more convenient to represent it in the form 
of constant energy loss rate (per unit time) as 

W\\ = / dii^ (13) 



n 00 

^ (zi,W||) = - j y dz2 j d qii (qii -vii) ^ 3me (zi, Z2; -q|| • V||) 



In the above, we used the odd parity of the imaginary part of the inverse dielectric function. The rate of the energy 
loss expressed per unit length is usually referred to as the stopping power S — . 

B. Charged Particle moving perpendicular to the surface (v y = 0) 

Now, let us assume that the particle has only velocity component perpendicular to the x^y plane, i.e., — 0. We 
shall focus on the last integration in Eq. ([8]), by changing the variable w = gzWz, so that we can perform the series 
string of simplifications on it, namely, 

1 r Qiu}{z2-z-i_)/v^ 

^—^ , , ,2 *^"^ (14) 

1 f g2"(z2-Zl)/-«z I r Q-^i^{z2-Zl)/v^ 

^ — duj uj— — — -26^^ {zi,Z2;q\\,u}) 2 , , ,2 i''*y^ (21, 2:2; QII, t^) 

" ' ql + {uj/v,y " ' vjj q^^+{oj/v,y 





00 

2i f duj UJ 



vtJ q^ + {uj/vz) 
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Consequently, by combining Eqs. ([14]) and ([8]), we obtain the following energy loss in the perpendicular case as 



2 oo oo oo 

K) = ^%^5m / dz, [ dz, f d\\ f e-' (zi, z^; 9|| , a;) e^^^-^^)/''^ . (15) 

(27r) eQ J J J J (qnv^) + 

— oo — oo ^ " ' 

Now, let us apply the energy loss formalism to single and double layers of of epitaxially grown graphene. The epitaxial 
form of graphene usually features the energy gap, which may be controlled by an external electric field from a gate. 



C. Charged Particle traversing/reflected from the graphene layers at an arbitrary angle 

For the case when the charged particle crosses the layers, none of the components of the velocity is zero. As in the 
perpendicular case, the initial position of the particle is of little importance and we may change variable zi + v^t — >■ zi 
in Eq. ([7]), followed by adding the obtained equation to Eq. ([8]), giving 



Wg{^)^~-^^^^ [ dzi [ dz2 [ dqj d^Ci\\e''i^^'--'^^e-^z,,Z2;q\\,Ci-^)q-\-w (16) 
(27r) w^en J J J J 



(27r) w^eo 
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2[Ze) 
(27r)^ v^tQ 



— OO — OO — OO 

OO OO OO 



3m j dzi J dz2 J dq, J d\\\e"''^'^-''h-^ {zi,Z2;q\i,(l-v)q~\-v 



-oo 



The last expression is obtained by changing the sign q — >■ — q and utilizing the symmetry relation Eq. (jlip . We now 
rewrite the energy loss via the energy- loss probability spectral function V (w). To do so, we introduce w = q • v = 
QzVz + q|| ■ V||, thereby obtaining 



We (v) ^ J W-P(a;,v) (17) 



2(^e)2 7 7 f exp(z(w-q|| •V|,)(z2-zi)M) 



oo oo 



2(ZeY f f f 

V{uj,v)^- — -3 — 5m / dzi / dz2 / d^q|| (zi, Z2; q|| , w) 

(2tt) €() J J J 



(27r)"eo J \J "J ^ ^' ' (q||^,,)V(cj-q|| .v||) 



—00 —00 



An alternative derivation of the above equation is provided in AppendixlAl One may show in a straightforward way 

that the above general expression may be expressed as Eq. (IT5|) in the case of V|| = 0. The other limiting case Eq. 

is less trivial. It is obtained from Eq. (|17p by taking the following limits lim — q|| • vy) /v^ = qz, lim dzi/vz = dti 

and using Eq. ^ . 

It is straightforward to derive the energy loss when the charged particle is elastically reflected from the Si surface 
(at zi = 0) of the SiC substrate, yielding 



00 00 



\ 4(Ze)^ f f /■ 2 _w cos((a;-q|| ■ {z2 - zi)/v,) 

V(uj,^r) ^ 5m / dzi / dza / ^ q|| e (zi, Z2; q|| , w) .2 , ^2"^ ■ (1^) 

(2^) eo J [qiiv,) + (cj - q,, • v,,) 

For multiple graphene layers separated by distance d, we can assume wave function localization on the layers, and 
neglect their overlap. Therefore, we obtain— 1^ 

e-i (zi,Z2;g||,o.)) =,5(zi-Z2) + ^i;,e'««l^i-^'^lnf (<5,- - nf i;,e-'ii l^-^'l'^) (z2 - /d) , (19) 

where the 2D Fourier transform of the Coulomb interaction has been introduced by Vq = ltsq\\ with = 47reoeb 
for background dielectric constant £(,• The arguments (g||,a;) of the noninteracting polarization function II^'^'' on the 
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j*^ layer have been omitted. Equation (flQl) can be substituted into Eq. (fT7|l and, after calculating the integrals over 
Zi and Z2, we obtain: 



Vr (w, v) 
Vn (w, v) 



4(Ze) 
(27r)^o 



I 7 



(^pz)^ + - q|| • V||)' 



r 
7^ 



(20) 



Additionally, we have introduced auxiliary expressions for reflective TZ and transmissive T particle trajectories: 



r = g|| cxp[-iq^(j - j')d] , 

n - 2<7i| cos [q,{3 - 3')d] + e^'ll^'* [q,sln{qj d) ~ q\\ COs(g,/d)] 
as well as the random-phase approximation (RPA) polarization matrix elements: 



(21) 
(22) 



For either a single (upper) or double (lower) layer configuration, Eq. ((20|) assumes the form 



(23) 



Pr(TC) (w, v) = 

(27r) eo 



(g||W,)' + (w-q|| -vii)^ 



1 - vqn\ 



(0) 



r(7^) [(i - v,nf^) (i - - vi e-^n^ 



Here, we have redefined the transmitted (refiected) notations as 



(24) 



r = 2q|| vl e-^i^" cos(g.d) n'^^nf + V, (n(°) + - 2v, uf^K 



7^ = e-39ii^w^ (n 



(0) «9||<1 



[2911 



^2g||d 



92 sin 



{2q^d) - g|| cos(2(2'2d)] 



+n(°) {«,nf [q, (1 ~ e^""'^) sin(<z,d) + <Z|| (Se^""'^ - l) cos(g,d) - 4q|| e 
+e29ii'^ [2(j|| e"""^ + s±n{q,d) - q\\ cos{q,d)] }) . 



(25) 
(26) 



III. ENERGY LOSS IN MLG 



Let us consider MLG with the SiC{0001) substrate lying in the z = plane. We shall ignore contribution to EELS 
from cr— electrons, substrate surface imperfections and its optical (FK-) phonons— as well as the resonant plasmon 
coupling with single-particle excitations.— For simplicity, we shall limit ourselves to only one or two graphene layers 
located sX z = d and z = 2d. In accordance with the classification in Ref.fl^ those are referred to as OML and IML 
correspondingly. One of the effects of the substrate is n— doping of the graphene layers, thus causing the chemical 
potential jii > 0, where i — 0,1 labels the layers. Hereafter we presume that both layers are held at the same chemical 
potential fi. The chemical potential /i is related to the Fermi wave vector kp — fi/hvp. The other effect is more 
subtle, and affects only zero layer. This layer, often referred to as a "buffer layer", exhibits crossover from Dirac 
to conventional 2DEG by virtue of the substrate induced gap Eg. Ultimately, this gap is related to the symmetry 
breaking between A and B sublattices of the buffer graphene layer. This effect depends on two main factors, namely, 
the angle between the T — K line and the principal axis of the SiC substrate and the degree of the hydrogen passivation 
of the substrate. Formally, the tt— electron dispersion around the K [K') point becomes 
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Ek = ^{hvpkf + {Eg/2f , (27) 

where Wi? is the Fermi velocity for free standing graphene. The effect of the substrate on the other graphene layers 
is mitigated by the buffer layer and the electrons obey conventional Dirac dispersion {Eg = in Eq. ([27|) ). Since it 
plays a crucial role in EELS (|24|) . we give the full form of noninteracting polarization along the real frequency axis: 

uf\q,io + iO+) = ^ + , (28) 

X {[iG>(a;i,_) - iG>(a;i,+)] 1< + [G<(a;i,_) + iG>(xi,+)] 2< 
+ [G<(xi,+) + G<(xi,_)] 3< + [G<(xi,_) - G<(a;i,+)] 4< + [G>(a;i,+) - G>(xi,_)] 1> 
+ [G>(a;i,+) + ^G<(xl,-)] 2> + [G>(a;i,+) - G>{-xx,-) - 1^2 - xl)] 3> 
+ [G>(-a;i,_) + G>(xi^+) - i7r(2 - x^)] 4> + [Go(xi,+ ) - Go(xi^-)] 5>} . 

Here, the following notations for the region functions have been introduced: 



Xo = 



El 



X2,± 



2fi± hu 



X3 = Jh^vj,q^+E^ 

G<(x) = 



G>(x) 



with the corresponding regions defined as 



Xq 


— x^ 


x^ 


Xq 


x^ 


Xg 



1< = 6 {fl — X2,- — hid) , 

2^—0 { — hw — /Lt + X2,-) 6* (/iW + /i — X2,-) {/J + X2,+ — ftw) , 

3< = 9 + X2^- - hu) , 

4< = 9 (hu + fi ~ X2.+) (hvpq — huj) , 

1> = 9 {2kF - q) 9 {huj - X3) 9 {fi + X2- - huj) , 

2y — 9 {hu — fi — X2^-) 9 {11 + X2,+ ~ huj) , 

3> = (/iw — /i — 2;2_-|-) 

4> = 9{q-2kF)9{haj-X3)9{fi + X2.--hoj) , 

5> = (fto; — hvpq) 9 (x3 — /icj) . 



Given Eq. ([28|) . the plasmon dispersion relation can be obtained as zeros of the real part of Eq. (jl9p . Alternatively, 
those solutions are given by the poles of the imaginary part of the RPA polarization in Eq. ([^5]) . The regions of (w, q) 
space of non-zero imaginary part of non-interacting polarization (|28p are referred to as particle-hole continuum. The 
plasmon branches in the particle-hole regions are Landau damped, and if excited rapidly lose their energy to the 
single-particle excitations. 

The plasmon branches for kpd = 1 and several values of Eg are presented in Fig.[TJ In the long-wavelength limit 
the single graphene layer (OML) exhibits a single undamped plasmon branch [Fig.[lja)] given by 



ujI = qV{Eg) 



(29) 
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FIG. 1: (Color online) Undamped (two higher frequency branches) and damped (two lower in frequency branches) plasmon 
dispersion relations for gapped graphene. The curves in (a) correspond to OML(single layer), the curves in (b) to lML(double 
layer) configuration. The red, green and blue curves show the plasmon dispersion for Eg/ = {0.0, 1.0, 1.5}, respectively. 



where the plasmon Drude factor is defined as 



V{E,) = ^^[l-{Ej2^,f] . (30) 

In the double layer configuration (IML), there are two plasmon branches, i.e., the symmetric w-i- and asymmetric lo- 
modes. For large interlayer distance kpd ^ 1, the two branches are qualitatively similar and given by 



w+ = qDiQ) , (31) 
Lol=qV{Eg) . 

In the opposite hmit of closely placed layers with kpd <C 1, the asymmetric branch becomes acoustic with frequency 
defined by 



ujI = 2dq^V{Eg) . (32) 

We note that the asymmetric branch is always smaller in frequency than symmetric mode. In the next section, we 
shall combine Eqs. (IMl) . (P5|) with (1^51) and simulate the energy loss for two distinct cases, namely when the charged 
particle motion is parallel and perpendicular to the graphene layers. 



IV. NUMERICAL RESULTS AND DISCUSSION 



We begin our discussion with the energy loss rates in Eq. (|T3| for single and symmetric double layer configurations. 
Our energy loss results for charged particle motion parallel to the surface are shown in Fig.[5J The separation between 
the two layers was chosen as kpd — 1.0. The charged particle travels above the graphene surface at a height kpzi — 2.0. 
All frequencies were analytically continued to the complex plane and acquire a small positive loss rate w — > w + 17 
with 7//i = 10~^. We notice that the plasmon contribution becomes smaller with increasing energy gap Eg. This 
can be easily attributed to decreasing plasmon Drude factor 'D{Eg). The contribution from the particle-hole modes 
is largely unaffected by the gap. In the double layer configuration, a small spike at the beginning of the plasmon 
contribution region is due to the two branches. Both undamped optical and acoustic plasmon branches contribute. 
However, with increased particle velocity, only the optical branch contributes, as it is demonstrated in Fig.[TJ Since 
the spectral weight of the optical branch is larger than that of the acoustic branch, we see that the spike appears. 
We note that one can obtain a closed-form analytic expression for the loss rate if we include only the long wavelength 
plasmon contributions.— In any event, this must be the leading contribution for small and moderate particle velocities 
v\\/vp < 6, because of the e~'^^^ /q factor. 

We now consider a single, free-standing (Eg ~ 0) and double-layer graphene and present results for the energy 
loss spectra of a charged particle moving perpendicular to the graphene surface. By symmetry, the transmission and 
refiection spectra from the single layer must be identical. The transmission spectra are shown in Fig.|31 The energy 
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loss was calculated using either the full version of the noninteracting polarization function in the RPA ([28)1 shown in 
Fig-El^a), or its long- wavelength plasmon pole approximation [Eqs. ([^ and ([?T|) ] as seen in Fig.[31Jb). One may classify 
three different scattering regimes depending on the charged particle velocity Vz/vp- In the low-velocity regime, as 
shown in column (1) of Fig.[3]for single and double layers, the particle-hole excitations dominate the scattering. There 
is an absorption spike due to the Landau damped plasmon mode which separates the particle-hole continuum from 
the undamped plasmon modes as shown in Fig.[TJa). The strength of the spike diminishes with increased charged 
particle speed. The position of such resonance absorption at (w 1.5) is independent of Vz, indicating the effect due to 
linear dispersion of the damped plasmon branch. In the intermediate velocity regime given in columns (2) and (3) of 
Fig-El the damped and undamped plasmon mode contributions are comparable. However, in the high-velocity limit, 
shown in columns (4) of Fig.|31 the undamped plasmons dominate the spectra. In that regime, one can employ the 
polarization in its long wave limit to a good approximation, as we can see in Fig.[31[b) . The approximation works for 
both the single and double-layer configuration. In the work of Allison, et al.,~ an attempt was made to compensate 
for the red shift of the absorption maximum by introducing a restoring force in the electron liquid. However, the 
origin of this restoring force was not clear. Additionally, the approximation introduced in their work does not yield 
results which converge on those when the RPA polarization for graphene is employed. Another observation is that 
the double layer just doubles the absorption of the single layer, without any perceptible spectral shift. This is a 
consequence of large interlayer distance. When the layers are closely packed the absorption peak is blue shifted with 
increasing number of layers, in agreement with Ref.l22l 

In Fig.m we present the transmission and reflection spectra for epitaxial MLG. The zeroth layer only acquires the 
energy gap. This gap increases in columns (1) through (3). Our principal observation is that in the intermediate 
velocity regime, given by columns (2) and (3) in Fig. 21 the absorption spectrum splits into two peaks. The one 
identified with the damped plasmon peak {huj/ii ~ 1.5) comes from the upper graphene layer for which Eg — 0. The 
other peak may be attributed to the symmetry-broken zeroth layer absorption (Eg ^ 0). We note that in the high 
velocity limit, shown in column (4) of Fig.|4l the gap results in a red shift of the absorption maximum. The small 
splitting of the peak is also visible on the damped plasmon contribution. As a matter of fact, the layers are so far 
apart that one cannot identify a plasmon mode as the acoustic branch. The reflection spectrum qualitatively mimics 
the transmitted spectrum, but doubles in height. This is because the charged particle spends twice as much time in 
between the layers compared to the transmissive case. Of course, the particle expends most of its energy on this part 
of the trajectory. 

To see the acoustic branch (Fig.[TJb)) in the spectra, we vary intralayer distance in (rows (a) through (c) of Fig.[S] 
for all velocity regimes [panels (l)-(4)]. Several values of the gap are also shown on the graph. The low-frequency 
acoustic branch boosts the long wavelength absorption (huj/fi w 0). This is especially so for small energy gap, since 
the separation between the branches is well pronounced [Fig. [Tib)]. In the high velocity regime or result os small 
interlayer separation [Fig. [5] (a. 4)] ore result is qualitatively similar to Fig. 3 of Ref.H^. With increasing angle of 
incidence Oi = taji~^(v\\/vz) the sharp peak of the particle hole absorption becomes less pronounced and finally 
reachesS^ the smooth curve as in Fig.[2][a.l- a. 4). The plasmon peak is blue shifted with sin^j finally reaching those 
of Fig. [21a. 5). Since the plasmon momentum is proportional to sin^^ one can interpret this blue shift as the Raman 
shift due to plasmon activation.'^*' Given its linear nature we deuce that for small interlayer distance mostly acoustic 
plasmon branch is activated, thus explaining Fig. 2(d) of Ref.ll^. This interpretation is further confirmed by the 
comparison with the large interlayer separation Raman shift, which follows V sin^i pattern. Unfortunately, as it 
follows from Fig.[T](b), it is hard to deduce existence of the gap from the acoustic plasmon branch. Experimentally 
one would have to transfer MLG onto a non-polar substrate and observe the increase in the plasmon slope if the gap 
was originally present. 

Comparing the plasmon absorption in the parallel and perpendicular cases, the plasmons contribute more to the 
former in the low- velocity regime, whereas in the latter case, their contribution is most pronounced for higher incoming 
charged particle velocities. 

V. SUMMARY 

We have investigated the role played by a gap in the energy dispersion on the absorption spectra of single and double 
layer configurations. All velocity regimes for the external charged particle moving perpendicular to the graphene 
surface were reported in Figs.[3l[4] and [H The plasmon pole approximation for the polarization agrees well with the 
results obtained with the full polarization in the RPA in the high-velocity regime only. The speed of the external 
charged particle determines whether the plasmon or particle-hole excitations dominate the scattering. Consequently, 
since gapped graphene has a different plasma excitation spectrum than free-standing graphene, its stopping power 
may carry distinct signatures of the substrate induced gap. We also demonstrated that our formalism can qualitatively 
describe experimental data of ELLS. It also allowed us to interpret the observed linear plasmon dispersion as coming 
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from the acoustical undamped branch. 
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Appendix A: An alternative derivation of the energy loss formula 

Here, we derive Eg. pT)) by a method similar to that in Refs.[25land[26l. First, we introduce the 2D Fourier transform 
and its inverse given by 



oo 



(Al) 



(A2) 



where we have adopted the notation of Camley and Mills with a — 1,2. By applying such a Fourier transformation 
defined in (jA2| to the Poisson equation (jSj, we obtain after some straightforward algebra the following differential 
equation 



g«("2-q2||-V|i)z2/fz 



Vz^O 



By applying the boundary conditions $ext (q2||; ±00, W2) — 0, the solution may be written as 

^ey^t (q2|h22,a;2) = - 



{l2\\Vzy + (t^2 - q2|| • V||)^ 



Therefore, the external potential assumes the form 



*ext {TC2\\,Z2,t2) ^ J^2d,2^ext (q2|| , ^2 , ^2) 



ZeVz 



d.2 e--*^ / d\2, -"^ ■ ^ ' " ■ ""^ 



{2t:Y Cq J " ('72||Wz) + (^2 - q2|| • V||) 

Similarly, we make a Fourier representation for the nonlocal inverse dielectric function as 

(ri||,^i,^i;i"2||,^2,i2) = (zi, 2:2; ri|| - r2||,^i - ^2) 

00 

(27r) J J 

— CO 

Combining Eqs. (IA5|) and (jA6|) with Eq. ([2]), we obtain 



(A3) 



(A4) 



(A5) 



(A6) 



*tot (ri||,zi,ti 



(27r)'eo 



duj e-"^'' I d^qii e*ii|-ri|i / dz: 



_ exp (i (w - q|| • V||) Z2/V, 

xe (zi,Z2;q||,w^ 



{q\\Vz 



' + (w-q|| -vii)^ 



(A7) 
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Here, we used the following identities: 



J d2r2|| e-'('5ii-i-ll)-"-=ii i2n f S (q,, - q2||) . 

The force acting on the particle is given by combining Eq. (jA7p and Eq. ([T]), with the replacement Vi = d/dvti. The 
resulting expression is inserted into Eq. (jH), yielding 



— OO — OO — OO — OO 

exp(i (w-q|| •V||)2:2M) 



ce 1 (zi,Z2;q||,w) 2 . s2 



00 00 



{Zef f f [2 [,. -If N exp[i(w-q|| •V|j)(2;2-2;i)M] 
— ^ / dzi / dz2 / d q|| / dwiwe (zi,Z2;q||,t^) . 2 , ^2 ^ 

(27r) eo^^ _i ^ _i ((zr.) +(^-qrv||) 



Now, let us change the dummy variables ui — > —lo qy — > — q|| in the above equation and utilize the symmetry relation 
in Eq. dn]). We obtain 



00 00 



^.(v)--^^^n. / dz, f dzJd\Jd..e-^[z,.z,,c,,,.)'-^^^ (A9) 
(2^) eo _J J J J (gp.) + (t^-q,| -vii) 
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in the sense of a; ~ v^^ sinOi. 

In our theory, we force the particle on the prescribed trajectory, thus making the incidence and the scattering angle to be 
the same. Therefore, Eq. ([T]) of Ref. 22 confirms that the plasmon momentum q is proportional to the sine of the incidence 
angle. 
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FIG. 2: (Color online) Stopping power, as a function of charged particle velocity in units of (^e/27r)^e7^ . Rows (a), (b) 
correspond to double and single layer configurations, respectively. Panels (1), (2) and (3) correspond to Eg/n = {0.0, 1.0, 1.5} 
illustrated by red, green and blue curves, respectively. The solid black curve is the particle-hole contribution, and the dashed 
curve shows the plasmon contribution. Row (4) gives the particle-hole contribution for the three values of the gap, and row (5) 
shows the associated plasmon contributions. 
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FIG. 3: (Color on-line) Energy loss spectra of the charged particle transmitted through single and double {kpd = 100) free- 
standing graphene. Rows (a) and (b) correspond to the full version of the RPA polarization and its long-wave plasmon 
approximation, respectively. In columns (1) through (4), the charged particle speed increases as v^/vf = {0.5,2.0,5.0,10.0}. 
The thin black curves correspond to particle-hole (and damped plasmon) contributions only. 
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FIG. 4: (Color on-line) Energy loss spectra of the charged particle transmitted through and reflected from a double (fefd ~ 100) 
epitaxial graphene. Rows (a) through (c) are for increasing energy gap on the zeroth layer with Eg/ ^ = {0.5, 1.0, 1.5}. In 
columns (1) through (4), the charged particle speed increases as v^/vf ~ {0.5, 2.0, 5.0, 10.0}. The thin black curves correspond 
to particle-hole (and damped plasmon) contributions only. The color schematic correspond to Fig. [1] and Fig. [2] 
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FIG. 5: (Color online) Energy loss spectra of the charged particle transmitted through a double free-standing graphene. Rows 
(a) through (c) correspond to increasing interlayer distance kpd — {5, 10, 50} . In columns (1) through (4), the charged particle 
speed increases as with v^/vp = {0.5, 2.0, 5.0, 10.0}. The thin black curves correspond to particle-hole (and damped plasmon) 
contributions only. Red thick and Blue dotted curves stand for Eg/ ji — 0.0 and Eg/ = 1.5 correspondingly. 



